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N omenclat ure 


A Hardening parameter 

[Bl] Matrix relating the incremental linear strain to the incremental displacement 

[-Bjv] Matrix relating the incremental non-linear strain to the incremental 

displacement 

[C] Elastic stiffness matrix 

[C] EP Elasto-plastic stiffness matrix 

C EP i Elasto-plastic stiffness tensor 

E Young’s modulus 

{F} Static force vector 

K hardening coefficient 

[K] Stiffness matrix 

[Ki] Linear part of the stiffness matrix 

[Af.v l] Non-linear part of the stiffness matrix 

Ni Shape function of the element 

[Q\ Matrix containing the shape functions of the element 
Sij Piola-Kirchoff stress tensor 

S t Surface of the body on which the stresses are specified 

S u Surface of the body on which the displacements are specified 

V Volume of the body 

{a} Flow vector 

{de} Incremental Green-Lagrange strain vector 

detj Incremental Green-Lagrange strain tensor 

dt Incremental time 

d u Incremental displacement 

{de} Incremental natural strain vector 

d tij Incremental natural strain tensor 

{da} Incremental Cauchy stress vector 

dcrij Incremental Cauchy stress tensor 

{/} Internal force vector 



Nomenclature 


Vll 


hi Initial thickness of strip 

h 2 Final thickness of strip 

n Hardening exponent 

r % reduction 

t Time 

u Displacement vector 

x Position vector 

A-rjij Incremental non-linear strain 
e§ 9 Equivalent plastic strain 
H Shear modulus 

v Poisson’s ratio 

ttij Angular velocity tensor 

4> n Neutral point angle 

p density 

{cr } Deviatoric stress vector 

ay Yield stress 

cr eq Equivalent stress 

<7ij Cauchy stress 

Deviatoric part of o ZJ 
<jij Jaumann stress rate 



Abstract 


To analyze the residual stresses in rolling, under plane strain condition, an updated 
Lagrangian elasto-plastic finite element formulation is presented. Inertial forces and 
roll deformation have not been considered. The material is assumed to yield according 
to the von-Mises criterion. Power-law type of strain hardening is used. Effect of 
temperature and strain rate on hardening behavior have been neglected. Friction at 
the roll-strip interface is incorporated only indirectly, as the boundary conditions at the 
interface are expressed in terms of incremental displacements. Parametric studies for 
change in reduction and the ratio of roll radius to strip thickness have been presented. 
The distribution of residual stresses is presented both over the domain (in the form of 
contours) as well as over the thickness. 



Chapter 1 


Introduction 


Rolling can be considered as one of the most important processes in industry, and is 
at the same time one of the oldest metal-working processes. Strip rolling is a process 
in which thickness of a flat sheet is reduced by passing it between two counter-rotating 
cylinders which are separated by a fixed distance. The process may be performed at 
room temperature on a cold mill or at higher temperature on a hot mill, both of which 
may consist of one or more mill stands. In cold rolling, the roll radius is usually more 
than 50 times the initial strip thickness. If the width of the strip is at least five times 
the length of the arc of contact, the non-plastic material prevents the lateral spread, 
and the deformation takes place effectively under plane strain condition. 

The complete analysis of rolling process is complex. As the strip enters the roll 
gap, first it deforms elastically. The relative velocity between the roll and the strip is 
such that frictional forces draw the metal in. The criterion of plastic flow governs the 
manner in which the transformation from elastic to plastic takes place. As the strip 
proceeds through the roll gap, more plastic flow occurs, until it comes out of the roll 
gap. The strip is then unloaded elastically. It is observed that the relative velocities of 
the roll and the strip change during the process. At the entry, the strip moves slower 
than the roll. As the strip starts accelerating forward, it reaches the surface velocity 
of the roll at the ” no-slip” or neutral point, where the roll and strip move at the same 
speed. As the further compression occurs, the strip speeds up and the direction of 
friction changes in such a way that it now retards the motion. Tensions are often 
applied in cold rolling in order to ensure adequate control of the dimensional accuracy 
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of the rolled strip. 

In rolling, the plastic deformation is often non-uniform. Therefore, to maintain the 
compatibility, different elements of the sheet cross-section unload to different levels of 
elastic deformation after they come out of the roll-gap. This induces self-equilibrating 
internal stresses in the sheet, called as residual stresses. It is generally believed that 
bad shape or poor flatness of rolled sheet is caused by residual stresses. Additionally, 
residual stresses have significant effects on further processing of the sheet. Thus, the 
distribution of residual stresses across the cross-section of cold rolled sheet plays an 
important role in rolling technology. The traditional method fail to predict residual 
stress distribution, because of the assumption of rigid-plastic material. The residual 
stresses can be taken into account by elasto-plastic finite-element analysis. 


1.1 Literature Survey 


To analyze the rolling process, numerous investigations, both analytical and experi- 
mental have been carried out. In view of the complexity of the process, analytical 
methods make use of several assumptions. These methods, therefore, provide useful 
but limited information. 

Finite Element technique has been applied to analyze the rolling process in recent 
years. Most studies have used the Eulerian formulation, as it is found to be convenient 
for steady-state processes like rolling. Further, it takes much less computational time 
compared to the other formulation, namely updated Lagrangian formulation. Both 
the rigid-plastic as well as elasto-plastic constitutive equations have been used. A rigid 
plastic analysis in the Eulerian reference frame is found to be most suitable when the 
primary concern is to find the roll force, torque and inter- facial stresses. However, 
elastic effects must be taken into account when the objective is to find the plastic 
boundaries and residual stresses. 

Eulerian elasto-plastic formulations have been reported by many authors. Dawson 
and Thompson (1978) have developed “initial stress rate” method. In this approach, 
the stress and velocity fields associated with the connective terms are assumed known 
and placed on the right-hand side of the finite element equations for momentum. An 
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iterative procedure is then. used during which these terms are corrected until there is 
no discrepancy between the new stress distribution and the previous one. Variations 
are performed on pressure and velocity only. Results on residual stress distribution or 
plastic boundaries have not been reported. 

Shimazaki and Thompson (1981) have used the mixed Eulerian formulation to solve 
for velocity, pressure and stress simultaneously. In their method, they have decoupled 
the momentum and continuity equations from the constitutive equation. The former 
set has been used to obtain velocity and pressure with stress assumed constant. The 
latter set has been used to obtain stress w r ith velocity and pressure held constant. 
Iterations between these two sets are carried out until convergence. Thompson and 
Berman (1984) have employed the procedure developed by Shimazaki and Thompson 
(1981) for the analysis of plane-strain rolling. For a perfectly plastic material, analysis 
has been done for cold rolling. However, rate of convergence has been slow and results 
presented are without getting full convergence. Effective stress contours are plotted, 
but individual stress components have not been presented. 

Abo-Elkhier et al. (1988) also have proposed an Eulerian formulation for solving 
elasto-plastic problem. For a linearly strain-hardening material, analysis for extrusion 
process has been conducted and some results on residual stresses are presented. 

Thompson and Yu (1990) have presented an Eulerian formulation to solve the 
rate equilibrium equation, in their steady state form, for velocity and stress. The 
authors have used the formulation for the analysis of strip drawing for a perfectly plastic 
material. It is noted that the rate of convergence deteriorates as the inelastic response 
becomes large. Effective stress contours are plotted without getting full convergence. 
No result is presented for residual stress distribution. 

Maniatty et al. (1991) have proposed an Eulerian elasto-plastic formulation in 
which the kinematics consists of a multiplicative decomposition of the deformation 
gradient into elastic and plastic components. Hart’s (1976) constitutive model which 
contains an internal scalar variable, has been used for solving rolling and extrusion 
problems. In rolling problems a simple sliding friction law, t t = P(v 0 - v t ) where 
ft is the traction on roll strip interface, v 0 and v t are the velocities of roll and strip 
respectively and /? is friction coefficient is used. The system of equations is highly ill- 
conditioned and the algorithm lacks stability as the material becomes less rate sensitive. 
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An improved algorithm for solving the problem has been presented by Maniatty (1994) 
which makes use of discontinuous pressure field and allows the unknown nodal pressure 
to be eliminated from the formulation. Only longitudinal stress variation for a typical 
case is presented in these references [Maniatty et al, 1991, ;Maniatty, 1994]. Results for 
other components of residual stresses and hydrostatic pressure are not presented. 

Malinowski and Lenard (1993) have proposed an elasto-plastic finite element scheme, 
using the Prandtl-Reuss constitutive relations, in an Eulerian-updated Lagrangian ref- 
erence frame. Problem is solved in two stages, first the velocity field is determined by 
minimizing the total power of deformation and then the stresses are obtained by path 
line integration of the Jaumann stress rate. Path line integration, though an efficient 
method, gives unrealistic results due to non-satisfaction of equilibrium condition. A 
less efficient approach which involves minimization of a functional under certain con- 
straints is used. Two types of shape functions (parabolic and Hermitian) are used 
and the results differ significantly from each other. In the proposed material model, 
Poisson’s ratio is assumed to be a continuous material function defined both for elastic 
and elasto-plastic behavior. Malinowski (1993) has used the methodology developed 
by Malinowski and Lenard (1993) for determining residual stress distribution in plane 
strain cold rolling, in which roll deformation has also been taken into account. An 
additional functional is minimized iteratively for the roll deformation. Only limited 
results are presented and computational difficulties are not discussed. 

Dixit and Dixit (1997) have reviewed various Eulerian elasto-plastic formulations 
and concluded that there are many unsolved difficulties in obtaining residual stresses 
using these formulations. They have proposed one simplified method to determine the 
more dominant longitudinal residual stress. 

As far as the updated Lagrangian formulation is concerned, there are some studies 
on rolling [Liu et al, 1985“ 1985 6 ] which use this formulation. But they have not 
studied the residual stress distribution. 
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1.2 Objective and Scope of the Present work 

The present work is an attempt to study residual stress distribution in a cold rolling 
process. A large deformation elasto-plastic finite element (FE) technique is used. Even 
though the Eulerian formulation is convenient for such analysis and takes much less 
computational time than the updated Lagrangian formulation, here, the latter formula- 
tion is used. This is because, there are many un-solved difficulties with the convergence 
of the Eulerian formulation, especially in the unloading zone. Jaumann stress rate and 
Green-Lagrange strain tensors are used respectively as the measures of objective stress 
and non-linear strain. Newton-Raphson technique is used to solve the non-linear in- 
cremental equations of the updated Lagrangian formulation. 

The material is assumed to yield according to von Mises criterion. Power-law type 
of strain-hardening relationship is used. The effects of strain rate and temperature on 
hardening behavior are ignored. Inertial forces, which have effect during the start up 
of the process, are not considered. Plane strain condition is assumed since the width 
to thickness ratio in strip rolling is usually more than five. Front and back tensions are 
not included in the analysis. For simplicity, rolls are considered to be rigid. Friction at 
the roll-strip interface is incorporated only indirectly, as boundary conditions on the 
interface are expressed in terms of incremental displacements. 

A large deformation elasto-plastic updated Lagrangian FE code is first validated 
by comparing the results on roll torque and roll pressure distribution with the corre- 
sponding results of Eulerian rigid-plastic formulation. Residual stress distribution is 
obtained both over the domain ( in the form of contours) as well as over the cross- 
section. Parametric study is carried out to study effects of reduction and the ratio of 
roll radius to strip thickness. 


1.3 Plan of the thesis 


Organization of the thesis is as follows: 

In the second chapter, the updated Lagrangian formulation with material and geo- 
metric non-linearities and the corresponding finite element formulation, are presented. 
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The boundary conditions and evaluation of the secondary variables (roll pressure, roll 
torque, residual stresses, equivalent plastic strain) are discussed, at the end of the 
chapter. 

In the third chapter, first the model is validated, then residual stress distribution 
is obtained and finally parametric study is carried out. 

In the fourth chapter, conclusions of the present work as well as suggestions for the 
further work are presented. 



Chapter 2 


Finite Element Formulation for 
Large Deformation Elasto Plastic 
Analysis 


In this chapter, the mathematical model of a large deformation elasto plastic problem 
is developed. The process is modeled as a steady state plane strain problem. The 
constitutive equation for the mechanical behavior of the work material is stated. The 
updated Lagrangian finite element formulation is developed. 


2.1 Updated Lagrangian Formulation 

In the study of deformation of a body subjected to external loading, often the original 
undeformed and unstressed configuration of the body is used for the formulation of its 
equation of motion. This is known as Lagrangian Formulation. This formulation is con- 
venient for small deformation which is the case in a majority of engineering problems. 
In such cases, the deformed configuration does not deviate much from the original one 
and hence the deformation can be described by an infinitesimal strain tensor, known 
as the engineering strain tensor, for which the strain-displacement relation is linear. 

On the other hand, for large deformation problems one has to use a finite strain 
measure which is expressed by a non-linear strain-displacement relation. Furthermore, 
the equations of motion, when expressed in the reference configuration, depend on the 
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deformation. Hence, for all large deformation problems, the Lagrangian Formulation 
proves to be cumbersome with the governing equations being difficult to solve. In 
such cases, one solves the problem using an incremental method known as Updated 
Lagrangian Formulation. In this formulation, it is assumed that the states of the stress 
and deformation of the body are known till the current configuration, say at time t. The 
main objective is then to determine the incremental deformation and stresses during 
an infinitesimal time step d t, i. e. from time t to t + d t. Here, the current configuration 
is used as the reference configuration for obtaining the incremental values. Unlike the 
Lagrangian Formulation, an incremental strain tensor is used. The methodology is 
particularly useful for elastic-plastic materials because the stress-strain relationship in 
such materials is usually expressed in an incremental fashion because of its dependence 
on the path of deformation. 


2.2 Governing Equations 

2.2.1 Strain-Displacement Relations 


Figure ( 2.1) shows a strip element in its current configuration (which is being used 
as the reference configuration) at time t and a possible deformed configuration at time 
t + A t. The vectors if and u t+At denote the displacements of a point P of the body 
at time t and t + At. respectively. The position vectors of P at the above times are 
depicted as and x tJrAt , respectively. For a very small value of At, one can represent 
it as dt. Thus , from the figure it can be seen that d u is the incremental displacement 
of P for the time step dt. For small incremental deformation, one can represent the 
deformation tensorially using the following linear incremental strain tensor: 


d €ij = - (d Uij + duj.i) 


(2.1) 


Note that de^ is not the incremental engineering strain tensor but the incremental 
logarithmic or natural strain tensor. 


If the incremental deformation is large, the incremental strain-displacement rela- 
tionship has to be non-linear, i. e. : 

dey = ~ (du, j -f- diij i -t- du^^du^j ) 


( 2 . 2 ) 
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The tensor de^ is known as Green Lagrange strain tensor. Henceforth, the incremen- 
tal displacement is assumed to be large and thus the non-linear strain-displacement 
relation given by eq. (2.2) is used. 


2.2.2 Constitutive Equations 

The plastic pan of the incremental stress-strain relationship is obtained from the plastic 
potential of the material. For a material, yielding according to the von Mises criterion 
in an isotropic fashion, the plastic potential is given by Owen and Hinton (1980). 


F(o ij: p) = Oegioij) - o y {p) (2.3) 

Note that 

F = 0 (2.4) 

represents the wield criterion. The plastic potential F depends on the Cauchy stress 
tensor <jy through it’s second invariant o eq called as equivalent stress and defined by 


.O / / \ I 

= hVijVij ) 2 


(2.5) 


where <jh is the ieviatoric part of a^. Further, F depends on the variable yield stress 
of the material. ~ y , through a hardening variable p. For the case of strain hardening. 
p is identified as the equivalent plastic strain ef 9 , and hence defined as : 



( 2 . 6 ) 


and 



(2.7) 


Here, defy is the plastic part of the incremental strain tensor de^- and the integration 
in eq.(2.6) is to be carried along the particle path. The dependence of o y on p ( or e v eq ) 
is normally approximated by a power-law type of relationship 



( 2 . 8 ) 



Finite Element Formulation for Large Deformation Elasto Plastic Analysis 


10 


Here, K is called the hardening coefficient, n is called as the hardening exponent and 
<7y is the yield stress at zero plastic strain. 


The incremental plastic strain (de?-) is obtained from the plastic potential using 

the following relation : 

dF 

defj = dX^~ (2.9) 

OGij 

where dA is a scalar. This equation is called as the flow rule. Differentiation of eq.(2.3) 
with respect to <jy gives 


dF _ 3 , 

dOij 2<Je 9 <Ty 

Using this and from eq.(2.9), one can determine dA as: 


( 2 . 10 ) 


dA = de* 9 


(2.11) 


Further, the hardening relationship and the yield condition (eq. 2.4) can be used to 
express dA as : 


dA = !^ 

H 


d<7 g 9 


( 2 . 12 ) 


where 


dcr 

deV 


(2.13) 


the slope of the hardening curve, can be obtained from eq. (2.8). Substitution of eqs. 
(2.10) and (2.12) in eq. (2.9) leads to the following constitutive equation 


p 3 do^ 
ij 2 H o, 


eq 




(2.14) 


This constitutive relationship between the deviatoric stress tensor and the incremen- 
tal plastic strain tensor is not really convenient for an updated Lagrangian formulation 
for which the incremental stress-strain relationship is needed. This can be obtained 
from eq. (2.14) as follows: 


j _ 3 Ojj dOeq 

ij 2 Ho eq do kl 


d &ki 


(2.15) 


Note that, from eqs. (2.3) and (2.10), we get 


da eq _ dF _ 3 / 

dcr k i do hi 2o eq akl ' 


(2.16) 




Finite Element Formulation for Large Deformation Elasto Plastic Analysis 


11 


Substitution of eq. (2.16) in eq. (2.15) leads to the following incremental plastic stress- 
strain relationship : 


<H- = 




dcr w . 


(2.17) 


The incremental elastic stress-strain relationship is given by 

d e\j = -^[-vdcrkkSij + (1 + ^)d <7y] (2.18) 

where defj is the elastic part of dejj , E is the Young’s modulus and v is the Poisson’s 
ratio. Adding the two relationships, we get 


d eij = d e\j + de?- 


+ — 


(i + ^) 


n ' ' i 


E 




d o 


Jd 


( 2 . 19 ) 


This is the incremental elastic-plastic stress-strain relationship needed in the updated 
Lagrangian formulation. However, it is the following inverse relationship [Chakrabarty 
(1987)] which is more useful : 


d &ij — Ejjkideki 


where 


Cg£ = 2/i 


fiikSjl + ~ 


V 


fiijfiki — 


9F°ij°ki 


2v' v '~ 2(3 im + H)o%\ 

and n is the shear modulus related to E and v by the relation: 

E 




2(1 + v) 


( 2 . 20 ) 

( 2 . 21 ) 

( 2 . 22 ) 


In eq (2.20), if the increment consists of only a pure rigid rotation, then d must 
be zero. A stress increment which satisfies this condition is called an objective stress 
increment. Thus, in the incremental elastic-plastic stress-strain relationship (eq. 2.20), 
the incremental stress tensor must be an objective stress tensor. The Cauchy stress 
tensor does not satisfy the objective requirement. So dcqj in eq. (2.20) can not be 
interpreted as the incremental Cauchy stress-tensor. The two most commonly used 
objective stress tensors are (i) the 2nd Piola-Kirchoff stress tensor 5^, (ii) stress tensor 
whose material time derivative is given by the Jaumann stress rate Cij. If we use the 
latter, then dmj in eq. (2.20) should be identified as df. Thus 

on dt = C^deH 


(2.23) 
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The Jaumann rate <j°j is related to Cauchy rate a\j by the relation 

o°ij d t = o'ijdt - a ik (Q jk dt) - o jk (Q lk dt) (2.24) 

where the incremental rotation tensor (Qydi) is defined as 

Dijdt = ^(dttij - d (2.25) 

Equations (2.21, 2.23-2.25) represent the elastic-plastic incremental stress-strain rela- 
tionship to be used in the updated Lagrangian formulation. 

2.2.3 Equilibrium Equations 

In a steady state problem, inertial forces are negligible. Neglecting body forces also, 
the equilibrium equation becomes 

<Tijj = 0- (2-26) 

This equation is to be satisfied in the deformed configuration (i. e. in the configuration 
at time t+dt) which is not known yet. In updated Lagrangian formulation, we need 
equilibrium equation in terms of the incremental stress. The incremental form of eq. 
(2.26) is given by Hill (1950). 


dttjjj (i7ij tk )du kt j — 0. (2.2 1) 

Note that duy in the above equation is the incremental Cauchy stress tensor and hence 
equal to <x°j dt. 

2.2.4 Incremental Governing Equations 

It is assumed that the strip passes through the roll gap in several increments, and 
after a certain number of increments, we arrive at the current configuration (i. e. the 
configuration at the present time i), which is treated as the reference configuration (see 
fig. 2.1). Further, the displacement vector u, and the stress tensor cry at the present 
time are known. Now a incremental displacement is applied to the boundary of the 
current configuration. Then to find the resulting incremental displacements, strains 
and stresses (dut, de^ , da tJ ), we have to solve the following equations: 
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(a) Incremental equilibrium equations: 

do'ijj = Oj 


(2.28) 


(b) Incremental stress-strain relations: 

Oij d t = C|£de w ; (2.29) 

[ < Tij defined by eqs. (2.24-2.25) and C^i defined by eq. (2.21) ] 

(c) Incremental strain-displacement relations: 

d eij = — (dtijj T du jtl -f- du^jdu^j ) (2.30) 

2.3 Finite Element Formulation 

In this section, a finite element formulation corresponding to the incremental equations 
of the updated Lagrangian approach of section 2.2 has been developed. First, an elastic- 
plastic matrix relating incremental stress and strain vectors is derived using the plastic 
potential proposed in section 2.2.2. Then, the equilibrium equations in the deformed 
configuration are derived using the principle of virtual work following a procedure given 
in Bathe et al (1975). Finally, the numerical scheme required to solve the finite element 
equations is discussed. 

2.3.1 Incremental Stress-Strain Relationship Representing 
Material Non-linearity 


The incremental stress-strain relation can be expressed in vector form as 



{<M = [C] EF { de} 

■ (2.31) 

where 


{d(j} d<J X y, d(7 Z z} 

(2.32) 

and 


{de} T = {de IX , de yy , 2de Iy , de z2 } 

(2.33) 
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are respectively the vector forms of the incremental (objective) stress and strain tensors. 
The expression for the elastic-plastic matrix [C EP ] is derived from the plastic potential 
F. For its derivation, it is convenient to express the yield condition eqs. (2.3) and (2.4) 
as 

F({o},p) = a eq ({a}) - a y (p) = 0 (2.34) 

where 


{o - } {&XXi &yy, &xyi Fzz} 


(2.35) 


On the yield surface, dF = 0 and hence 


or 


, dF Tf " dF J 

(_) {<M + = 0 


{a} r {dcj} — ^4dA = 0 


where, the flow vector {a} is defined by 


w 


dF 

5{a}’ 


the parameter A is given by 


, 1 dF, 

A = ~d> ;^ dp ' 


(2.36) 

(2.37) 

(2.38) 

(2.39) 


and dA is the same scalar which appears in eq. (2.9). Eq. (2.10) gives the following 
expression for the flow vector: 

{a} = (2.40) 

^ eq 

where {a}' is the vector form of the deviatoric part of Oij. Eqs. (2.6) and (2.11) give 

dp = dA = de%. (2.41) 


Substitution of eqs. (2.34) and (2.41) in eq. (2.39) leads to the following expression for 
A: 

H/t_. r\rr 

(2.42) 


day = d^ 
dp deeq 


The total strain increment can be split into elastic and plastic parts. Thus 

dF 


{de} = {de e } + {de p } = [C\- x {do) + dA 


d{c} 
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= [C]" 1 {da} + dA{a} (2.43) 

Here, [C\ is the matrix form of the elasticity tensor Cijki (See eq. 2.18). The second 


part of the right hand of eq. (2.43) is due to the flow rule (eq. 2.9). 
both sides of equation (2.43) by {a} r [C'] and using eq. (2.37), we get 

Premultiplying 


{a} T [C\{de} 

A + {°} r [C]{a} 

(2.44) 

Substituting the above expression for dA in eq. (2.43), we get 



,, ■> (v-ri-irj i { a } T [C]{de} , . 

{de} - [C] {da} + A + { G pp]{ G } W 

(2.45) 

which leads to 

r j -i |Vl [C]WW T [C] \ r . , 
d<7 “ ( lC] A + {a} r [C]{a} j <d6} 

(2.46) 

Therefore, 

, n sp (, n [CK“}{“} T [C] \ 

lC1 -(> C1 A + {anC]{a}) 

(2.47) 


where the expression for {a} is given by eq. (2.40). For a material with power strain 
hardening law (eq. 2.8), the expression (2.42) for A reduces to 


A = Kn{e%) n ~ l 


(2.48) 


For an isotropic material, the expression for [C] for plane strain case is given by 


[Cl 


(1 + u)(l - 2u) 


' l-v 

V 

0 

U 


V 

1 - V 

0 

V 

( 2 - 

0 

0 

1-2 v 

2 

0 

V 

V 

0 

1 - V 



As stated at the end of subsection 2.2.2, the vector {da} in eq. (2.46) should be 
interpreted as an objective stress increment. 

2.3.2 Updated Lagrangian Formulation for Geometric Non- 
Linearity 

When a body experiences a large deformation and/or rotation, the equilibrium must 
be established in the current configuration. As stated earlier, the normal solution 
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procedure adopted in this case is updated Lagrangian. That is to say that when all 
the kinematical variables are known from time 0 to time t in discrete time steps, the 
objective is to establish the equilibrium in the configuration at time t + At. In order to 
derive the finite element equilibrium equations at time t + At , the principle of virtual 
work requires that 

J V t + At a £ At dT/i+At = Rt+At (2.50) 

where crf+ At is the Cauchy stress tensor at time t + At , 5(e)£ At ) is the virtual infinites- 
imal strain tensor at time t + At , R t+At is the virtual work of the surface forces at 
time t + At (body forces being neglected) and V t+At is the domain at time t + At . 
The major difficulty in application of eq. (2.50) is that the configuration at t + At is 
unknown. Moreover quantities like Cauchy stress tensor are not purely additive as one 
has to take care of rotation into account. 


An elegant way of formulating the problem is given by Bathe et al (1975). Here, the 
virtual work expression (eq. 2.50) is transformed to an integral over the known domain 
at time t (i. e. V t ) using the 2nd Piola- Kirchoff stress tensor and the Green-Lagrange 
strain tensor. Thus, the virtual work expression becomes 

J t S$ Ai 6{ t e% At ) dH* = R t+At . (2.51) 

where Sy denotes the 2nd Piola- Kirchoff stress tensor whose relation with the Cauchy 
stress tensor <jy is given by 

(2.52) 

[The right superscript stands for the current configuration and the left subscript denotes 
the reference configuration.] 


Since the Cauchy stress tensor is always referred to the current configuration, the 
left subscript has been omitted. Here t+ A tA,m represents the derivative of the position 
vector x 1 at time t with respect to the one ( x t+At ) at time t + At and denotes 
the ratio of densities at time t and t + At. The Green-Lagrange strain tensor t e(£ At is 


defined as 


teU At 


_ £ ( t + At , t+A t , t+At t+At'l 
— 9 \ tU io + tU 3,i + tU k,i t U k,j ) 


(2.53) 


where tu\j At is the derivative of the displacement vector u t+At at time t + At with 
respect to the position vector at time t. 
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An incremental decomposition of stress and strain gives 

t s*; At = tSfj + t ASij 




(2.54) 


St€y At = 5 t Atij + 5 t Ar)ij , 

(2.55) 

where, 

t^Cij = — (t A + 1 A Uj,i) 

(2.56) 

and 




“ (^A Uk,i) Au^j) . 

(2.57) 


Here t^ u i,j stands for the derivative of t A u ( incremental displacement vector at time 
t ) with respect to x* (position vector at time t ). Thus, eq. (2.51) with incremental 
decomposition can be written as 


[ tASySitAeiijdV* + [ tASySUArHjdV* + [ 4 « 5 ( t Ar^dV 4 

J'V * J\At 

+ J cl/it Ae^dV 1 = R t+At . (2.58) 

The 2nd integral on the left hand side is a higher order term compared to other terms 
and hence can be neglected. Moreover, using an elastic-plastic incremental stress- 
strain relationship (eq. 2.20), tASij can be approximated as tC^i tAtki ■ Here, the 
left subscript of tC^i denotes the time at which it is to be evaluated. The above 
simplification will definitely give rise to an error in the right hand side of eq. (2.58). 
The error is normally neutralized by using some iterative scheme like Newton Raphson 
Method. Now, eq. (2.58) can be written as 

,Cf£, A£ H i(,Aey)dK‘ + l v / ij &(Ani 3 )&V t + 

= R t+At . (2.59) 


Owing to the symmetries in tAe^, t A and ajj, eq. (2.59) can be rewritten 
in the following form: 

J vt (a,{A<} T ) t {C] EF t { Ae)dV‘ + J vt (i5 i {At 7} t ) [T}\{ A,}dV* 

+ J (5,{Ae} T ) {a}‘dV‘ = R ,+At . 


(2.60) 
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where, 

i{Ae} = t Ae yy , 2 t Ae xy , t Ae zz } T , (2.61) 

{a} 1 = {cr‘ xx , o l yy , a l xy , <4} T » ( 2 - 62 ) 

{tAu XiX , jAj-~.j. tAu XjZ , tA u ytX , tAu y , y , tA u y , z , } , (2.63) 

[r ; = 

and 

Ef- 

The t[C)“ p matrix is the elastic-plastic matrix evaluated at time t which is given by 
eq. (2.47 . 

2.3.3 Finite Element Equations 

The domain is discretized into a number of elements and the incremental displacement 
field over a typical element e is approximated as 

t{An} = [Q] t {A uY (2.66) 


!EP 0 
o E]‘ 


(7 * 

xx xy 

°yx a yy 
0 0 


0 

0 


(2.64) 


(2.65) 


where, 

t{Au} e = {, At- . t A«J, . . . , t At£, * A u*}. 

and t Aul . tAu y stand for the (unknown) incremental displacements of node V in x 
and y directions respectively . Here 


where 



(2.67) 


{Qi} T = [ Ak 0 N 2 0 ... N n o], 

{Q 2 } r =[o .V, o N 2 ... 0 N n ]. (2.68) 


and Ni, the functions of ( x, y ), are called the shape functions. 



Finite Element Formulation for Large Deformation Elasto Plastic Analysis 


19 


With displacement field defined by eq. (2.66), the strain field within the element 
can be calculated in terms of nodal displacements. From eq. (2.66), 

^ = t{Q,} T j «{Au}* (2.69) 

where 

< m 5 = dJ § f - (2-70) 

Substituting eq. (2.69) in eqs. (2.56) and (2.57) and arranging them in vector form, we 
get the strain displacement relations: 

*{Ae} = t [B L ] *{Au} e , ( 2 . 71 ) 


t{ At?} = t [B N ] t {A< 


where, 


t[B L ) = 


t{Qi}J x 

<{Q 2 } T y 

tWS+aoi} 

0 


and 


( [Bn] T = [ t{Ql},x t{Ql},y 0 t{Q 2 },x i{Q 2 },y 0 j 


(2.72) 


(2.73) 


(2.74) 


Using the strain displacement relations derived above, the contribution to the inte- 
gral (eq. 2.60) from a typical element e (with volume V* ) can be written as 

MM' T ( J vt t[Bi] T «[Cf p ,\B l ] dV ) ,{A .«}• + 

(j,, i[B«f PT ABn] «") + (2.75) 

<5,{Au}' r (7 ,[B l ] t {<t}‘ dV*) = M Au}‘ t (,+a,{F}‘) 

Here, the contribution to the term R t+At from the element e has been expressed in terms 
of the elemental external force vector t+A t {F} e using a standard procedure. Since the 
variation in the incremental displacement vector <U{ Au} e is arbitrary, expressing the 
term within 1st parenthesis as linear elemental stiffness [ Kj_ ] e and that of 2nd one as 
nonlinear elemental stiffness [K NL ] e the above equation can be written as 


(t[Ki] e + t[i^wi] e ) i{A.-u} e + t{f} e — t+Af-fF 1 } 6 


(2.76) 



Finite Element Formulation for Large Deformation Elasto Plastic Analysis 


20 


or 

t[I<] e t{Au} e + t {fY = t+ At{F} e . (2.77) 

where 

«{/}' = ! ,\Bl] T {ct>‘ dV>. (2.78) 

J Kt 

is the internal force vector. 

Assembling the elemental matrices t[K] e and the elemental vectors t{f} e and £ _ 
over all the elements, we get the following global equation : 

,[K] ,{A«} + ,{/} = 1+il {.F}. (2.79) 

Here £ [A] is the global stiffness matrix, £ {Au} is the global incremental displacement 
vector (at time t ) and £ {/} and t+±t{F} are the global internal and external force 
vectors respectively. Decomposing t +At{F}, eq. (2.79) can be written as 

t[K) £ {Au} + t {f} = t {F) + t {AF}. (2.80) 

Here, £ {F} is the (global) external force vector at time t and £ {AF} is the (global) 
incremental force vector from time t to t + A t. In updated Lagrangian formulation, it 
is assumed that the equilibrium equations have been satisfied exactly at time t. Thus, 

= (2.81) 

Then, eq. i 2.80) reduces to 

£ [A] £ {Au} = £ {AF}. (2.82) 

Now, suppose we find the incremental stresses corresponding to the solution ( £ {Au}) 
of eq. (2.82). add those to the earlier stresses to find the stresses at time t + At and 
use eq. (2.78) to determine the internal force vector £+ a t{f} at time t + At. it is 
possible that it may not match with the applied (or external ) force vector. Thus, 
eq. (2.82) represents only an approximate equilibrium equation at time t + At (The 
approximation is mostly due to the linearization and simplification involved in the steps 
between eqs. 2.58 to 2.59). - A solution of such an approximate equation may involve 
a significant amount of error and, depending on the time/incremental displacement 
step, may becomes unstable. Therefore, it is necessary to modify eq. (2.82) to turn 
it into an iterative problem capable of providing a solution with desirable accuracy. 
Newton Raphson scheme is one of the most commonly used technique to handle such 
nonlinearity. 
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2.3.4 Numerical Scheme 

Normally we encounter two types of problems in mechanics: (1) Load control problem 
and (2) Displacement control problem. In a load control problem, the desired deforma- 
tion is achieved by prescribing a load at a point. On the other hand, in a displacement 
control problem, it is the prescribed displacement which gives the desired deformation. 
In a load control problem, if the load falls due to change in geometry, then most iterative 
schemes fail to converge. Displacement control problems axe free from this limitation. 
In this thesis, only displacement control problems are considered. 

(I) The Newton Raphson Scheme 

The scheme is best described with the help of Fig. (2.2) in the following steps. In 
Fig. (2.2), m is an equilibrium point at load level {P} . The aim is to reach another 
equilibrium point ( m 4- 1 ) at a load level ( {P} + A{P} ). 

Step 1 — First of all, the stiffness matrix t [K] in eq. (2.82) is evaluated using the 
stresses corresponding to point m. Next, eq. (2.82) is solved with t{AF} = {AP} to 
get an incremental displacement vector ( {Au}. 

Step 2 — From the incremental nodal displacements, incremental strain components 
are calculated by using eq. (2.2). Incremental stresses are calculated from the incre- 
mental strains by using suitable constitutive relation. Generally incremental stresses 
are calculated at the Gauss points of the elements. If the point under consideration is 
in elastic range, the elastic stress-strain matrix [C] is used. Otherwise, elastic-plastic 
matrix [C) EP is used. Since elastic-plastic matrix is dependent on the current state 
of stress, the stress value at position m is used for the computation of elastic-plastic 
matrix. Eq. (2.23) is used to calculate Jaumann stress rate and the total Cauchy stress 
follows from eq. (2.24). Step 3 - Once the stresses are calculated at all the Gauss 
points of all the elements, the internal nodal force vector t +At{f} is calculated. This 
force corresponds to the equilibrium load at point i in Fig. (2.2). 

Step 4 - Since the applied nodal forces ( t {P } + t {AP} are not equal to internal 
nodal forces t +A t{f}, the deformed configuration is not an equilibrium configuration at 
load level ( t {P} + t {AP}). It gives rise to an unbalanced nodal forces (t{P} + t{AP} - 

t+Ai{/})- 
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Step 5 — The steps 1 to 3 are repeated with unbalanced nodal forces as the load 
vector and i as a starting point. Successive repetition of the above steps will lead 
to the equilibrium point (m + 1) for load level t {P} + t {AP} when the unbalanced 
nodal forces become almost zero(less than a preassigned quantity). In each successive 
iteration .incremental displacement vector is added to its previous value. 

In full Newton-Raphson method, the stiffness matrix is updated after every itera- 
tion. This takes a lot of computational time per iteration but the number of iterations 
required for convergence, in one increment, is quite small. When one tries to reduce the 
computational time per iteration, it invariably increases the number of iterations per 
increment required for convergence. This may actually result in increasing the overall 
computational time per increment. The best compromise between the reduction in 
computational time per iteration and the increase in number of iterations per incre- 
ment is offered by the modified Newton-Raphson method in which the stiffness matrix 
is modified only after the increment is over. The iterative scheme of the modified 
Newton-Raphson method can be described by the following equations: 


t [K] t {Au}U = t {P} (i - 1} (2.83) 

for i=l,2,3.... 

where t [K) is the same as in equation (2.82) and the right side vector is given by 

,{B} (01 = <+a.{A F) (2.84) 

= 1+at {AF} - , +a ,{AF}«> (2.85) 

for i=l,2,3.... 

However, this scheme creates some convergence problems for the rolling process. In 
this process, incremental deformation is simulated by prescribing incremental transla- 
tion of the strip. While doing so, it must be ensured, in every increment, that a node 
exists at the entry to the roll gap. Thus, the minimum distance through which the strip 
can be advanced in one increment is the distance between two consecutive nodes on 
the top face. To achieve the convergence using the modified Newton-Raphson scheme, 
the incremental deformation has to be sufficiently small, but this will require a large 
number of elements and consequently a large number of dof. Solving a problem of such 



Finite Element Formulation for Large Deformation Elasto Plastic Analysis 


23 


a large size on IIT K computing systems takes a lot of time. To keep the problem at 
a manageable level, the element size has to be sufficiently big. This means the incre- 
mental deformation becomes large leading to convergence problem. To overcome this 
difficulty the following strategy is adopted. In this scheme, t [K] in equation (2.83) is 
replaced by the stiffness matrix updated after carrying out the first iteration. Then 
the modified Newton-Raphson method is applied using this new stiffness matrix. 

(II) Scheme for divergence 

If, during the iterations of modified Newton Raphson scheme, the unbalanced force 
vector in (i + 1)^ iteration increases compared to the unbalance force vector of i t h 
iteration, it means the scheme is diverging. This happens because we get a large value 
of t Au^ I+1 ^. In this situation, if magnitude of the stiffness matrix is increased and 
(i + l) th iteration is repeated, we can get a smaller value of t Au^ +1 \ and convergence 
can be achieved. If again there is divergence, the magnitude of the stiffness matrix can 
be increased again. 

One way to increase the magnitude of the stiffness matrix is to add to it a fraction 
of the elastic stiffness matrix [K] E corresponding to that configuration. It is calculated 
using only the linear part of t [K] and replacing [C] EP by [C] E . Thus, the modified 
stiffness matrix becomes 


t[K} new = (l-a) t [K] + a[K) E (2.86) 

In this work, the above scheme is used to arrest the divergence, where a is found 
by numerical experiment. 

(Ill) Euler Forward Integration Scheme In step 2 of modified Newton-Raphson 
scheme, if equation (2.31) is used to calculate the incremental stress vector where [C KP ] 
is based on the current equilibrium configuration, then a lot of error occurs in estimating 
the stresses leading to a large value of unbalance. For better estimation of stresses, 
usually Euler Forward integration scheme is used. This scheme is described below. 

• If the state at time t is elastic, 

^ strain increment using eq. (2.2) and Jaumann stress increment assum- 
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ing elastic behavior. 

‘{Act} = [C*] t {Ae} (2.87) 

2. Evaluate Cauchy stress increment using eq. (2.24). Calculate t+At {cr} 

3. Determine t c eq and t+At a eq using eq. (2.5) 

4. If t+At (J e q < t+At &y then elastic behavior holds. [If t+At a eq — t+At o y , the Gauss 
point is neutrally loaded]. RETURN. 


5. If 
culate 


t~\- At - 


' eq 


> 


t+At 


<j y a transition from elastic to plastic state has occurred. Cal- 


t+At r 


Ratio = 


J eq 


t/X 

u eq u eq 


( 2 . 88 ) 


and change the state to plastic. Calculate the Cauchy stress using the following pro- 
cedure. 


In each iteration, the sub-incrementation method is followed and the stress tensor 
components are updated after each sub-increment by adding the sub-incremental stress 
components corresponding to the elasto-plastic sub-incremental strain. First the C EP j 
corresponding to the last updated state is calculated. Next, Jaumann stress increment 
for the sub-increment is evaluated. Then Cauchy stress increment is evaluated using 
eq. (2.24). Finally, the updated state of stress is evaluated: For i th sub-increment, 
the Jaumann increment is given by 

adt (i) = t+At [C EP } {i - l) d t { Ae} (2.89) 

and the updated state of stress is obtained as 

t+At {o} {i) = f+At {< 7 } (i-1) + {da} [i) (2.90) 


In the above equation, 

t+At {a}° = t {a} + Rqtio[C E ] t { Ae} 

t+At^EPjO = t^pEP^ 

and d t {Ae} is obtained from the sub-incremental displacement 


d t {Au} = 


(1 - Ratio) t {Au} 


n 


(2.91) 

(2.92) 

(2.93) 
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where n is the number of sub-increments. 

Use above equations to find t+At Oij. RETURN. 


• If the state at time t is plastic, the sub-incrementation method in (5) above is 
applied with Ratio set to zero. 

2.3.5 Unloading 

The unloading has to be incorporated if one wants to get residual stresses in elasto- 
plastic deformation. The unloading criterion defined by Chakrabarty (1987) is imple- 
mented. According to this criterion,. unloading occurs if 


{n} T {Aa} < tolerance (2.94) 

where {Aa } is incremental stress vector and {n} represents the unit outward normal to 
the yield surface at the current stress point t {a} in a 9 — D space. Neutral loading and 
positive loading are represented by replacing the "less than” symbol in above equation 
by the ” equal to” and ’’greater than” symbols respectively. This unit normal to the 
yield surface is defined by 

{n} = d\{de p } (2.95) 

where dX (> 0) is given by equation (2.11) and {de p } is incremental plastic strain 
vector. 

When unloading is detected, the stress tag of Gauss point is changed from plastic 
to elastic and the increment is repeated. It is possible that during the repetition of the 
increment, the same Gauss point may again become plastic, but in that case it is not 
checked for unloading. 

2.3.6 Choice of Shape Functions and Numerical Integration 

Since the finite element method is a numerical procedure involving discretization of 
the domain, important considerations pertain to the convergence of the numerical 
solution with respect to mesh refinement. For monotonic convergence, the elements 
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must be complete and the elements and mesh must be compatible. If these conditions are 
fulfilled, the accuracy of the solution increases continuously as we continue to refine 
the finite element mesh. This mesh refinement should be performed by subdividing 
a previously used element into suitable number of more elements; thus, the old mesh 
must be ” embedded ” in the new mesh. The requirement of completeness of an element 
means that the polynomial approximating the displacements of the element must be 
able to represent the rigid body displacements and the constant strain fields. The 
requirement of compatibility means that the displacements within the element and 
across the element boundaries must be continuous. 

Eq. (2.59) suggests that the convergence criterion will be satisfied if d u x and du y 
are chosen to be bilinear functions of x and y. To satisfy the above condition, we 
have chosen the second order Serendipity approximation with eight nodded rectangu- 
lar elements for the finite element formulation. Here, the term x 2 y 2 of the assumed 
approximation polynomial is omitted because it increases the size and the bandwidth 
of the stiffness matrix without improving the approximation outside the element. The 
approximation is normally expressed in terms of the natural coordinates f and p. Thus, 
the shape functions TV) in eq. (2.68) are the second order Serendipity shape functions 
in f and p. 

Exact evaluation of integrals appearing in element coefficient matrices and right side 
vectors is not always possible because of the algebraic complexity of the coefficients in 
differential equations. In such cases, it is natural to seek numerical evaluation of these 
integral expressions. Numerical evaluation of integrals, called numerical integration, 
involves approximation of the integrand by a polynomial of sufficient degree, because 
the integral of a polynomial can be evaluated exactly. 

The most used numerical integration method is Gauss-Legendre numerical integra- 
tion and formula for integrals defined over a rectangular master element SIr is given 
by 


[ F(€,p) dp = f [ f F(£,p) dp df 

m n 

(2.96) 

« EE F(€i,pj) Wi Wj 
i=l i=l 

(2.97) 


where m and n denote the number of quadrature points (or the Gauss points) in the 



Finite Element Formulation for Large Deformation Elasto Plastic Analysis 


27 


£ and r] directions, (£$, ijj) denote their natural coordinates, and Wi and Wj denote the 
corresponding weights. 


2.4 Boundary Conditions 


The typical boundary conditions can be classified as: 

(i) Geometric or essential boundary conditions. 

The incremental displacement vector dui is specified at a point on the boundary. 

(ii) Force or natural boundary conditions. 

The incremental stress vector dt* = da^ nj is specified at a point on the boundary, 
where nj represents the components of the unit outward normal vector. 

Sometimes, we have mixed boundary conditions, i. e. at a point on the boundary, 
some components of incremental displacement vector and the remaining components 
of the incremental stress vector are specified. 

The specific boundary conditions used in the rolling problem (Fig. 2.3) are given 
below. 


1. The top free surface before entry to roll gap (AB): As stated earlier, 
the ro llin g process is simulated by specifying incremental displacements of 
some points on the strip. The points chosen are the 2 nodes on the boundary 
AB preceding the node B. At these points the boundary condition are 


d u x — du* x 

Essential 

(2.98) 

# P 

II 

X) 

Essential 

(2.99) 


for node m, du* is the difference in x coordinates of node B and node m and 
du* = 0. For node n, du* x and du*. are calculated from the condition that 
it should be on the roll surface. 

Since AB is a free surface, at the remaining points of AB, the stresses can 
be specified to be zero. It is observed that no bulging occurs except near the 
roll entry. Therefore, to avoid bulging, the vertical incremental displacement 
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must be specified to be zero. Therefore, the following two conditions are 
chosen. 

d^ = d<j X y = 0 Natural (2.100) 

d u y = 0 Essential (2.101) 

2. The roll-strip interface ( BC ): The nodes on this boundary travel along 
the roll surface. Therefore the following boundary conditions are used. 


d u x = du* Essential (2.102) 

d u y = du* Essential (2.103) 

where du* x and du* are calculated from the mesh geometry. 

3. The top free surface after exit from roll gap(CD): The boundary CD 
is a free surface. So, the stress components can be specified to be zero, but 
in order to avoid bulging of the boundary, y component of the incremental 
displacement is specified as zero. Thus, the two boundary conditions are 


d t x = do xy = 0 Natural (2.104) 

d u y = 0 Essential (2.105) 

4. The front boundary ( DE ): The boundary DE is a stress-free boundary. 
Therefore boundary conditions are 


o 

II 

H 

b 

xJ 

III 
4* 

T3 

Natural 

(2.106) 

d ty — ^ 

Natural 

(2.107) 


5. The axis of symmetry, z-axis (EF) : Because of symmetry, the y- 
component of incremental displacement and ^-component of incremental 
stress vector are zero on EF. Therefore, the boundary conditions on EF 
are 

dt x = dcr xy = 0 Natural (2.108) 

d u y = 0 Essential (2.109) 

6. The Entry Boundary (FA): Boundary FA is a stress free boundary. 
Therefore, the boundary conditions are 


dtc — do xx — 0 

dty — d O x y 0 


Natural 


( 2 . 110 ) 

( 2 . 111 ) 


Natural 
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2.5 Evaluation of secondary variables 

As explained above, solution of the rolling problem is obtained in the form of nodal 
displacements and Gauss point stresses by solving equation (2.83). From these two sets 
of primary quantities, the secondary quantities like roll pressure, roll torque, residual 
stress distribution and contours of equivalent plastic strain are obtained residual stress, 
the shear and, normal stress distribution at the roll-strip interface, the stress and strain 
contours and roll torque are calculated. 

The details of these calculations are explained briefly in the following lines. 

i. Roll pressure or inter-facial normal stress: 


Using the values of Oij at the Gauss points, first they are extrapolated to the nodes 
at the roll-strip interface, then normal stress t n is evaluated as 



t n — tiUi 

(2.112) 

where 

~t% — GijTij 

(2.113) 

and rij is unit outward normal to the interface. 


ii. Roll torque: 



Roll torque is calculated as: 

T = AW/A9 

(2.114) 


where A 9 is incremental angle of rotation of roll and ATT is total work involved in 
the increment, which consists of the work of deformation and the work to overcome 
friction (front and back tensions assumed to be zero). 


Work of deformation is calculated using the Cauchy stress and linear incremental 
strain. Thus, 

. AW deformation = J y OyA£ijdV (2-115) 

Work to overcome friction is given by 

A W friction = [ t s (RA6 - Au s )ds 
J s 


(2.116) 
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where R is the roll radius, A u s is the tangential displacement and t s is the shear stress 
at the roll strip interface. The shear stress is calculated from 


l*»l = /*l*n| (2.117) 

where fj, is the coefficient of friction. The sign of t s is determined by using the following 
expression (with out the tensions) for the neutral point which is based on the slab 
method [Bland and Ford(1948)]. 

taU ' l {l 2 (i)n = ~ ^ ln {^} (2-118) 

where 

(j> n — location of neutral point. 
h\ — initial thickness of the strip. 
h\ = final thickness of the strip. 
r = percentage reduction. 


The parameter a is given by, 



(2.119) 


iii. Residual stress distribution: 

Distribution of residual stress components is first obtained over the domain in 
the form of contours. The contours are plotted from their Gauss point vaues using 
a standard technique. Then, the variation of residual stresses over the thickness is 
obtained from these contours. 

iii. Contours of equivalent plastic strain: 

These contours are also obtained from the Gauss point values using a standard 
technique. 
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Figure 2.2. Illustration of Newton Raphson iterative scheme for a single load 
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Chapter 3 


Results And Discussion 


The finite element modeling of cold strip rolling process developed in the previous 
chapter has been applied to various sets of input variables. The metal considered is 
hot rolled and annealed steel AISI 4140, the material properties of which are given in 
Table 3.1 [Reddy et al (1995)]. 


Young’ modulus : 210 GPa 

Poisson’s ratio : 0.3 

Initial yield stress : 418.9 MPa 

Hardening coefficient : 699.0 MPa 

Hardening exponent : 0.266 


Table 3.1: Material properties of AISI 4140 steel. 


The value of 0.08 is used for the coefficient of friction p, to determine the friction 
work in equation ( 2.116). 

First the model is validated. Then residual stresses and plastic strain are presented 
for one typical case in the next section. Parametric study is done by changing the 
parameters R/hi and % reduction. In all the cases, the initial thickness of the strip is 
kept constant. 

The mesh data for the various cases is given in Table 3.2 below. 
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reduction(%) 

R/h 1 

elements 

No. of nodes 

5 

130 

272 

961 

5 

65 

248 

877 

15 

65 

348 

1227 


Table 3.2: Mesh data for various cases. 


3.1 Validation of the model 


For validation of the model, the results are compared with the predictions of the Rigid- 
plastic Eulerian formulation of Dixit and Dixit (1996). Table 3.3 shows the comparison 
of roll torque with that found by Dixit (1997) based on the formulation of Dixit and 
Dixit (1996). It can be seen that the agreement is good with the error being 6%. 
It is to be noted that whereas the Eulerian formulation gives only the steady-state 
value of torque, the updated Lagrangian formulation gives the information about how 
the torque varies as the strip advances. This variation is exhibited in section 3.4 on 
parametric study (see Fig 3.12). 

Present work : 1.391 kN-m/m 

Dixit(1997) : 1.479 kN-m/m 

Table 3.3: Comparison of roll torque per unit width (r = 5%,i?/hi = 65). 


Figure 3.1 shows the comparison of normal stress distribution with that of Dixit 
(1997). The predicted values are lower than that of Dixit (1997) because the roll 
deformation has been neglected in the present work. Whereas the distribution is more 
or less uniform in case of Dixit (1997), the present work shows 2 peaks, the entry peak 
being smaller than the exit peak. This difference is also due to the assumption of rigid 
rolls. 


3.2 Residual stresses and strain 

Figures 3.2 - 3.4 show contours of stress components o xx & a xy and equivalent plastic 
strain e p q for a typical case of r = 5% and R/h\ = 65. It can be seen that, when 
the strip travels a sufficient distance (> 3 / 12 ) from the roll exit, the stresses and strain 
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remain more or less constant, except of course at the end because of the free surface 
effect. These represent the residual stresses and strain in the strip. 

Figures 3.5 - 3.7 show the variations of steady-state values of a xx k o xy and 
e eq along the thickness. It can be seen from Fig 3.5 that the longitudinal residual 
stress (&xx) is positive at the top surface whereas it is negative at the center. Similar 
trend has been reported by Dixit and Dixit (1997) who have used a simplified method 
based on Eulerian formulation to determine the residual stresses. Our values cannot be 
compared directly with their values as their results are for higher reduction. Figure 3.5 
shows that shear residual stress (c xx ) is very small compared to longitudinal residual 
stress (<r xx ). The equivalent plastic strain (e p eq ) is observed to be more at top fibers than 
at the center (Fig 3.7 ). This shows that the assumption of homogeneous deformation, 
often made in elementary method of analysis is not valid at least for the case of low 
reduction. 


3.3 Deformed Configuration and Plastic Boundaries 


Figure ( 3.8 a) shows the deformation pattern at the front end while Figure ( 3.8 b) 
shows the same at the roll exit. These figures show the deformed position of an initially 
straight cross-section of the strip. The pattern is different at the front end because of 
the free surface effect. The pattern at the roll exit is consistent with the distribution 
of longitudinal stress (cr xx ) along the thickness (Fig 3.5). Figure ( 3.8 b) shows that 
the elongation of top fibers is less than that of central fibers. It means, the unloading 
has to be more for the central fibers to maintain the compatibility in steady-state. 
Whereas the unloading for the top fibers results in retaining a part of the positive 
elastic strain (e xx ), the unloading for the central fibers continues beyond the zero value 
so as to acquire the negative elastic strain (e xx ). This phenomenon explains ■why the 
residual stress are positive in the top fibers and negative in the central fibers. 

Figure ( 3.9) shows the plastic boundaries. The front boundary is the first nonzero 
contour of e p and indicates the onset of plastic deformation. The second boundary 

eq 

indicates the beginning of unloading. Beyond this boundary, there is no change in the 
values of e* and therefore the contour pattern remains the same. 

e</ 
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3.4 Parametric Study 

In this section, the effects of process variables of % reduction and R/hi ratio are 
analyzed on residual stress and strain, roll torque and normal stress at the roll-strip 
interface. While analyzing these effects, the initial thickness of the strip and the friction 
coefficient are kept constant. Figures 3.10 - 3.13 show the effects of % reduction and 
R/hi on a xx , enroll torque and normal stress respectively. 

3.4.1 Effect of Change in R/hi Ratio 

Figure 3.10 shows that even though the distribution pattern along the thickness re- 
mains the same, the values of longitudinal residual stress (a xx ) decrease with R/h\ ratio. 
This happens because the deformation (especially the change of length in longitudinal 
fibers) becomes more homogeneous with the increase in roll radius. The distribution of 
e p eq ( both the pattern as well as values), however, seems to be in sensitive with R/hi 
ratio (Fig 3.11). The roll torque increases with R/h\ ratio (Fig 3.12) because the 
arc-length increases with roll radius. Figure 3.13 shows that the pattern of normal 
stress distribution doesn’t change with R/hi ratio. However, the values decrease. This 
may be because, now, the deformation is more homogeneous. 

3.4.2 Effect of Change in % Reduction 

More reduction means more deformation. As a result, even though the distribution 
pattern remains the same, the values of o xx (Fig 3.10), e p eq (Fig 3.11), roll torque (Fig 
3.12) and normal stress (Fig 3.13) increase. 
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Figure 3.1: Comparison of normal stress distribution along 
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Figure 3.3: Contours of for one typical case of 5% reduction and R/hi - 65 
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of a*! along thickness, for one typical case of 5% reduction and R/hi — 65 
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Figure 3.9: Plastic boundaries 
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Figure 3.10: variation of o xx along thickness for various cases 
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Chapter 4 


Conclusions and Suggestion for 
Further Work 


4.1 Conclusions 


Present work is an attempt towards determination of residual stresses in plane strain 
cold rolling. Even though an Eulerian formulation is a natural formulation for a steady- 
state process like rolling, there are unsolved difficulties with the convergence of elasto- 
plastic Eulerian formulation. Therefore an updated Lagrangian formulation is used. A 
large deformation elasto- plastic finite element code has been developed for the deter- 
mination of residual stresses. The code predicts that the shear component of residual 
stresses (a xy ) is negligible which is along expected lines. The thickness wise distri- 
bution of longitudinal residual stress (<J XX ) matches qualitatively with the published 
results. Thus, the code is successful in evaluating a xx and a xy components of the resid- 
ual stresses. However, very high values of a yy are predicted which is contrary to the 
expectation. In fact cr w also should be of the order of a xy . The reason for high values 
of a yy is improper modeling of the boundary condition on the free surface after the 
roll exit. The condition used in the present work is du y = 0 whereas more appropriate 
condition is ty = 0. However, when this condition was tried, certain numerical diffi- 
culties were encountered, which could not be overcome during the available time. But, 
certain suggestions to overcome the problem are offered in the next section. 

JENTRAL UBRAR> 
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4.2 Suggestions for Further Work 

• Evaluation of a yy : 

As stated above, for proper evaluation of a yy , the boundary condition on the top 
surface (after the roll exit) should be ty = 0 instead of du y = 0. However, when 
this condition was applied, the vertical incremental displacements of the nodes 
in the vicinity of the roll become so large that they penetrated the roll. This is 
not possible physically. For such nodes, the increment has to be repeated using 
the condition on vertical incremental displacement which now would be equal to 
the distance of the node from the roll surface. The finite element code should be 
modified to incorporate this feature for the nodes leaving the roll gap. 

• Inclusion of Inertial forces : 

Inertial forces cannot be neglected during start up of the process. During the 
start-up, some non-zero stresses are developed in the entry region. Further, the 
roll torque during start-up is more than the steady-state value. These features 
can be simulated if inertia forces are included. Provision can be made to include 
the inertia forces only during start up, i.e. till the roll torque reaches a steady 
value. After that inertia forces need not to be included. 

• Friction boundary condition : 

In the present work, the friction coefficient has been used indirectly in the calcu- 
lation of roll torque. This is because the shear stresses at the roll-strip interface 
are not predicted accurately by the present model. For this purpose, the constant 
friction law can be used as one of the boundary condition. 

• Roll deformation : 

The prediction of normal stress distribution at the roll-strip interface can be 
improved by incorporating roll deformation. This can be done using Hitchcock s 
formula as it is simple. 
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